Phase separation of a driven granular gas in annular geometry 
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This work investigates phase separation of a monodisperse gas of inelastically colliding hard 
disks confined in a two-dimensional annulus, the inner circle of which represents a "thermal wall" . 
When described by granular hydrodynamic equations, the basic steady state of this system is an 
azimuthally symmetric state of increased particle density at the exterior circle of the annulus. When 
the inelastic energy loss is sufficiently large, hydrodynamics predicts spontaneous symmetry breaking 
of the annular state, analogous to the van der Waals-like phase separation phenomenon previously 
found in a driven granular gas in rectangular geometry. At a fixed aspect ratio of the annulus, the 
phase separation involves a "spinodal interval" of particle area fractions, where the gas has negative 
compressibility in the azimuthal direction. The heat conduction in the azimuthal direction tends to 
suppress the instability, as corroborated by a marginal stability analysis of the basic steady state with 
respect to small perturbations. To test and complement our theoretical predictions we performed 
event-driven molecular dynamics (MD) simulations of this system. We clearly identify the transition 
to phase separated states in the MD simulations, despite large fluctuations present, by measuring 
the probability distribution of the amplitude of the fundamental Fourier mode of the azimuthal 
spectrum of the particle density. We find that the instability region, predicted from hydrodynamics, 
is always located within the phase separation region observed in the MD simulations. This implies 
the presence of a binodal (coexistence) region, where the annular state is metastable. The phase 
separation persists when the driving and elastic walls are interchanged, and also when the elastic 
wall is replaced by weakly inelastic one. 

PACS numbers: 45.70.Qj 



I. INTRODUCTION 

Flows of granular materials are ubiquitous in nature 
and technology [lj]. Examples are numerous and range 
from Saturn's rings to powder processing. Being dissipa- 
tive and therefore intrinsically far from thermal equilib- 
rium, granular flows exhibit a plethora of pattern forming 
instabilities [2|, |3| . In spite of a surge of recent interest 
in granular flows, their quantitative modeling remains 
challenging, and the pattern forming instabilities provide 
sensitive tests to the models. This work focuses on the 
simple model of rapid granular flows, also referred to as 
granular gases: large assemblies of inelastically colliding 
hard spheres 0, [H, H, 0, H, Q ■ In the simplest version 
of this model the only dissipative effect taken into ac- 
count is a reduction in the relative normal velocity of 
the two colliding particles, modeled by the coefficient of 
normal restitution, see below. Under some additional as- 
sumptions a hydrodynamic description of granular gases 
becomes possible. The Molecular Chaos assumption al- 
lows for a description in terms of the Boltzmann or En- 
skog equations, properly generalized to account for the 
inelasticity of particle collisions, followed by a systematic 
derivation of hydrodynamic equations [l(J, EH E2] • For 
inhomogeneous (and/or unsteady) flows hydrodynamics 
demands scale separation: the mean free path of the par- 
ticles (respectively, the mean time between two consecu- 
tive collisions) must be much less than any characteristic 



length (respectively, time) scale that the hydrodynamic 
theory attempts to describe. The implications of these 
conditions can be usually seen only a posteriori, after 
the hydrodynamic problem in question is solved, and the 
hydrodynamic length/time scales are determined. We 
will restrict ourselves in this work to nearly elastic colli- 
sions and moderate gas densities where, based on previ- 
ous studies, hydrodynamics is expected to be an accurate 
leading order theory [J, [|| 0, H, [1| . These assumptions 
allow for a detailed quantitative study (and prediction) 
of a variety of pattern formation phenomena in granular 
gases. One of these phenomena is the phase separation 
instability, first predicted in Ref. | I 3| and further inves- 
tigated in Refs. [U, [H, [H, [T3, Eallll • This instability 
arises already in a very simple, indeed prototypical set- 
ting: a monodisperse granular gas at zero gravity con- 
fined in a rectangular box, one of the walls of which is 
a "thermal" wall. The basic state of this system is the 
stripe state. In the hydrodynamic language it represents 
a laterally uniform stripe of increased particle density at 
the wall opposite to the driving wall. The stripe state was 
observed in experiment [201 ] , and this and similar settings 
have served for testing the validity of quantitative model- 
ing [2l|, [22|, [H| . It turns out that (i) within a "spinodal" 
interval of area fractions and (ii) if the system is suffi- 
ciently wide in the lateral direction, the stripe state is 
unstable with respect to small density perturbations in 
the lateral direction [ID, [T^, [l6| . Within a broader "bin- 
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odal" (or coexistence) interval the stripe state is stable 
to small perturbations, but unstable to sufficiently large 
ones [3 [l^]. In both cases the stripe gives way, usu- 
ally via a coarsening process, to coexistence of dense and 
dilute regions of the granulate (granular "droplets" and 
"bubbles" ) along the wall opposite to the driving wall 
[l7l Il3 ] . This far-from-equilibrium phase separation 
phenomenon is strikingly similar to a gas-liquid transi- 
tion as described by the classical van der Waals model, 
except for large fluctuations observed in a broa d reg ion of 
aspect ratios around the instability threshold [18]. The 
large fluctuations have not yet received a theoretical ex- 
planation. 

This work addresses a phase separation process in a 
different geometry. We will deal here with an assembly 
of hard disks at zero gravity, colliding inelastically inside 
a two-dimensional annulus. The interior wall of the an- 
nulus drives the granulate into a non-equilibrium steady 
state with a (hydrodynamically) zero mean flow. Particle 
collisions with the exterior wall are assumed elastic. The 
basic steady state of this system, as predicted by hydro- 
dynamics, is the annular state: an azimuthally symmet- 
ric state of increased particle density at the exterior wall. 
The phase separation instability manifests itself here in 
the appearance of dense clusters with broken azimuthal 
symmetry along the exterior wall. Our main objectives 
are to characterize the instability and compute the phase 
diagram by using granular hydrodynamics (or, more pre- 
cisely, granular hydrostatics, see below) and event driven 
molecular dynamics simulations. By focusing on the an- 
nular geometry, we hope to motivate experimental stud- 
ies of the granular phase separation which may be ad- 
vantageous in this geometry. The annular setting avoids 
lateral side walls (with an unnecessary/unaccounted for 
energy loss of the particles). Furthermore, driving can 
be implemented here by a rapid rotation of the (slightly 
eccentric and possibly rough) interior circle. 

We organized the paper as follows. Section [IT] deals 
with a hydrodynamic description of the annular state of 
the gas. As we will be dealing only with states with a zero 
mean flow, we will call the respective equations hydro- 
static. A marginal stability analysis predicts a sponta- 
neous symmetry breaking of the annular state. We com- 
pute the marginal stability curves and compare them to 
the borders of the spinodal (negative compressibility) in- 
terval of the system. In Section HTTl we report event-driven 
molecular dynamics (MD) simulations of this system and 
compare the simulation results with the hydrostatic the- 
ory. In Section IIVI we discuss some modifications of the 
model, while Section [V] contains a summary of our re- 
sults. 



II. PARTICLES IN AN ANNULUS AND 
GRANULAR HYDROSTATICS 

The density equation. Let N hard disks of diameter d 
and mass m = 1 move, at zero gravity, inside an annulus 



of aspect ratio f2 = i? e xt/-Rint, where i? xt is the exterior 
radius and R- mt is the interior one. The disks undergo 
inelastic collisions with a constant coefficient of normal 
restitution ji. For simplicity, we neglect the rotational 
degree of freedom of the particles. The (driving) interior 
wall is modeled by a thermal wall kept at temperature 
To, whereas particle collisions with the exterior wall are 
considered elastic. The energy transferred from the ther- 
mal wall to the granulate dissipates in the particle in- 
elastic collisions, and we assume that the system reaches 
a (non-equilibrium) steady state with a zero mean flow. 
We restrict ourselves to the nearly elastic limit by as- 
suming a restitution coefficient close to, but less than, 
unity: 1 — fi <C 1. This allows us to safely use granular 
hydrodynamics [f| . For a zero mean flow steady state the 
continuity equation is obeyed trivially, while the momen- 
tum and energy equations yield two hydrostatic relations: 



V • q(r) +1 = 0, p = const , 



(1) 



where q is the local heat flux, / is the energy loss term due 
to inelastic collisions, and P = P(n, T) is the gas pressure 
that depends on the number density n(r) and granular 
temperature T(r). We adopt the classical Fourier rela- 
tion for the heat flux q(r) = — «VT(r) (where K is the 
thermal conductivity), omitting a density gradient term. 
In the dilute limit this term was derived in Ref. [ill ]. 
It can be neglected in the nearly elastic limit which is 
assumed throughout this paper. 

The momentum and energy balance equations read 



V • [«VT(r)] = I, p= const , 



(2) 



To get a closed formulation, we need constitutive rela- 
tions for p(n,T), n(n,T) and J(n, T). We will employ 
the widely used semi-empiric transport coefficients de- 
rived by Jenkins and Richman [24| for moderate densi- 
ties: 
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and the equation of state first proposed by Carnahan and 
Starling [25| 



p = nT(l + 2G) 



(4) 



where G = v{l - ff)/(l - v) 2 and v = n (ttgP/4) is the 
solid fraction. Let us rescale the radial coordinate by 
i?i nt and introduce the rescaled inverse density Z(r, 0) = 
n c /n(r,8), where n c = 2/ (-\/3eP) is the hexagonal close 
packing density. The rescaled radial coordinate r now 
changes between 1 and f2 = i? e xt/^?int, the aspect ratio 
of the annulus. As in the previous work [lfij], Eqs. 
([4]) and ([3]) can be transformed into a single equation for 
the inverse density Z(r): 



V ■ [F(Z)VZ] = AQ{Z) . 



(5) 
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The dimensionless parameter A = (27r/3)(l— fi) (Rint/d) 
is the hydrodynamic inelastic loss parameter. The 
boundary conditions for Eq. (© are 

dZ(l,9)/d9 = and V„Z(Q,6) = 0, (7) 

The first of these follows from the constancy of the tem- 
perature at the (thermal) interior wall which, in view of 
the constancy of the pressure in a steady state, becomes 
constancy of the density. The second condition demands 
a zero normal component of the heat flux at the elastic 
wall. Finally, working with a fixed number of particles, 
we demand the normalization condition 



Z-\r, 0)rdrd9 = irf(fl 2 - 1) , 



(8) 



where 
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is the area fraction of the grains in the annulus. Equa- 
tions ©-([8]) determine all possible steady state density 
profiles, governed by three dimensionless parameters: /, 
A, and f2. 

Annular state. The simplest solution of the density 
equation ^ is azimuthally symmetric (^-independent): 
Z — z(r). Henceforth we refer to this basic state of the 
system as the annular state. It is determined by the 
following equations: 

[rT{z)z'\ = rAQ(z), z'(fi) = 0, and 

r« (9) 

z^rdr = (Q 2 - l)//2, 

where the primes denote r-derivatives. In order to solve 
the second order equation © numerically, one can pre- 
scribe the inverse density at the elastic wall zq = z(£l). 
Combined with the no-flux condition at r = 17, this con- 
dition define a Cauchy problem for z(r) [l6l . hj. Solv- 
ing the Cauchy problem, one can compute the respective 
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FIG. 1: The normalized density profiles obtained from the 
MD simulations (the dots) and hydrostatics (the line) for il = 
2, A = 81.09, and / = 0.356 (equivalently, za = 2.351). The 
simulations were carried out with N = 1250 particles, [i — 
0.92, and _Ri nt = 22.0. Also shown is a typical snapshot of the 
syst 
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FIG. 2: The main graph: the marginal stability curves k = 
k(f) (where k is an integer) for Q = 1.5 and A = 10 4 (circles), 
A = 5 x 10 4 (squares), and A = 10 5 (triangles). For a given 
A the annular state is stable above the respective curve and 
unstable below it, as indicated for A = 10 4 . As A increases the 
marginal stability interval shrinks. The inset: the dependence 
of k m ax on A 1//2 . The straight line shows that, at large A, 
k nc A 1 / 2 



value of / from the normalization condition in Eq. ([5]). 
At fixed A and Q, there is a one-to-one relation between 
zq and /. Therefore, an alternative parameterization of 
the annular state is given by the scaled numbers zq, A, 
and f2. The same is true for the marginal stability anal- 
ysis performed in the next subsection. 

Figure [T] depicts an example of annular state that we 
found numerically. One can see that the gas density in- 
creases with the radial coordinate, as expected from the 
temperature decrease via inelastic losses, combined with 
the constancy of the pressure throughout the system. 
The hydrodynamic density profile agrees well with the 
one found in our MD simulations, sec below. 

Phase separation. Mathematically, phase separation 



manifests itself in the existence of additional solutions 
to Eqs. ©-© in some region of the parameter space /, 
A, and fi. These additional solutions are not azimuthally 
symmetric. Solving Eqs. I©-© for fully two-dimensional 
solutions is not easy [l3|. One class of such solutions, 
however, bifurcate continuously from the annular state, 
so they can be found by linearizing Eq. ©, as in rect- 
angular geometry [13l llq] . In the framework of a time- 
dependent hydrodynamic formulation, this analysis cor- 
responds to a marginal stability analysis which involves a 
small perturbation to the annular state. For a single az- 
imuthal mode ~ sin(k6) (where k is integer) we can write 
Z(r,9) — z{r) +e H(r)sin(fc#), where H(r) is a smooth 
function, and £< la small parameter. Substituting this 
into Eq. © and linearizing the resulting equation yields 
a fc-dependent second order differential equation for the 
function T[r) = T[Z(r)\ H(r): 

1 /fc 2 AQ'(Z) \ 

T * + r Tk -\y + ^(zT) k ~ m 

This equation is complemented by the boundary condi- 
tions 

r(i) = o and r'(fi) = o. (11) 

For fixed values of the scaled parameters /, A, and fi, 
Eqs. (Unj) and (jTTJ) determine a linear eigenvalue prob- 
lem for fc. Solving this eigenvalue problem numeri- 
cally, one obtains the marginal stability hypersurface 
k = k(f, A, fi). For fixed A and fi, we obtain a marginal 
stability curve k — k(f). Examples of such curves, for a 
fixed fi and three different A are shown in Fig. [2] Each 
k = k(f) curve has a maximum k max , so that a density 
modulation with the azimuthal wavenumber larger than 
kmax is stable. As expected, the instability interval is 
the largest for the fundamental mode k — 1. The inset 
in Fig. [2] shows the dependence of k max on A 1 / 2 . The 
straight line shows that, at large A, k max cx A 1 / 2 , as in 
rectangular geometry fig] . 

Two-dimensional projections of the (/, A, fi)-phase 
diagram at three different fi are shown in Fig. [3] for the 
fundamental mode. The annular state is unstable in the 
region bounded by the marginal stability curve, and sta- 
ble elsewhere. Therefore, the marginal stability analysis 
predicts loss of stability of the annular state within a fi- 
nite interval of /, that is at f min (A, ft) < f < f max (-A, fi). 

The physical mechanism of this phase separation insta- 
bility is the negative compressibility of the granular gas 
in the azimuthal direction, caused by the inelastic energy 
loss. To clarify this point, let us compute the pressure of 
the annular state, given by Eq. ©. First we introduce 
a rescaled pressure P = p/{n c To) and, in view of the 
pressure constancy in the annular state, compute it at 
the thermal wall, where T — Tq is prescribed and z(l) is 
known from our numerical solution for the annular state. 
We obtain 
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FIG. 3: Two-dimensional projections on the [A, / (SI 2 — l)/2] 
plane of the phase diagram at = 1.5 (the solid line), fi = 3 
(the dotted line), and fi = 5 (the dashed line). The inset 
shows more clearly, for fi = 3, that the marginal stability 
curve (the solid line) lies within the negative compressibility 
region (bounded by the dashed line). 



The spinodal (negative compressibility) region is deter- 
mined by the necessary condition for the instability: 
(dP/df) A n < 0, whereas the borders of the spinodal 
region are defined by (dP/df) A n = 0. Typical P(f) 
curves for a fixed fi and several different A are shown 
in Fig. |U One can see that, at sufficiently large A, the 
rescaled pressure P goes down with an increase of / at 
an interval fi < / < /a< That is, the effective compress- 
ibility of the gas with respect to a redistribution of the 
material in the azimuthal direction is negative on this 
interval of area fractions. By joining the spinodal points 
fi and (separately) fi at different A, we can draw the 
spinodal line for a fixed fi. As A goes down, the spin- 
odal interval shrinks and eventually becomes a point at 
a critical point (P c , f c ), or (A c , f c ) (where all the critical 
values are fi-dependent). For A < A c P(f) monotoni- 
cally increases and there is no instability. 

What is the relation between the spinodal interval 
(A, Ji) and the marginal stability interval (/ mm , / mQX )? 
These intervals would coincide were the azimuthal wave- 
length of the perturbation infinite (or, equivalently, k — > 
0), so that the azimuthal heat conduction would van- 
ish. Of course, this is not possible in annular geometry, 
where k > 1. As a result, the negative compressibility 
interval must include in itself the marginal stability in- 
terval (fmin, fmax)- This is what our calculations indeed 
show, see the inset of Fig. |3J That is, a negative com- 
pressibility is necessary, but not sufficient, for instability, 
similarly to what was found in rectangular geometry [la ] . 

Importantly, the instability region of the parameter 
space is by no means not the whole region the region 
where phase separation is expected to occur. Indeed, 
in analogy to what happens in rectangular geometry 
[3 Gil i phase separation is also expected in a binodal (or 
coexistence) region of the area fractions, where the annu- 
lar state is stable to small perturbations, but unstable to 
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FIG. 4: The scaled steady state granular pressure P versus the 
grain area fraction / for Q = 1.5 and A = 1.1 x 10 3 (the dotted 
line), A = 1.5 x 10 3 (the dash-dotted line), A = 3.5 x 10 3 
(dashed line), and A = 5 x 10 4 (the solid line). The inset 
shows a zoom-in for A = 3.5 x 10 3 . The borders /i and f 2 
of the spinodal interval are determined from the condition 
(dP/df) An — 0. The thick solid line encloses the spinodal 
balloon where the effective azimuthal compressibility of the 
gas is negative. 

sufficiently large ones. The whole region of phase sepa- 
ration should be larger than the instability region, and it 
should of course include the instability region. Though 
we did not attempt to determine the binodal region of 
the system from the hydrostatic equations (this task has 
not been accomplished yet even for rectangular geometry, 
except in the close vicinity of the critical point (HI), we 
determined the binodal region from our MD simulations 
reported in the next section. 

III. MD SIMULATIONS 

Method. We performed a series of event-driven MD 
simulations of this system using an algorithm described 
by Poschel and Schwager (26|. Simulations involved N 
hard disks of diameter d — 1 and mass m = 1. After 
each collision of particle i with particle j, their relative 
velocity is updated according to 

%j = Vij - (1 + n) {Vij ■ fij)hj , (12) 

where Vij is the precollisional relative velocity, and = 
rij/\fij\ is a unit vector connecting the centers of the 
two particles. Particle collisions with the exterior wall 
r = i?ext are assumed elastic. The interior wall is kept at 
constant temperature To that we set to unity. This is im- 
plemented as follows. When a particle collides with the 
wall it forgets its velocity and picks up a new one from a 
proper Maxwellian distribution with temperature To, see 
e.g. Ref. pages 173-177, for detail. The time scale is 
therefore dim/Ty) 1 / 2 = 1. The initial condition is a uni- 
form distribution of non-overlapping particles inside the 
annular box. Their initial velocities are taken randomly 
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FIG. 5: Typical steady state snapshots for N = 1250 and 
n = 6 (a); N = 5267 and fi = 3 (b), and N = 6320 and 
Q = 6 (c). Panels (a) and (b) correspond to annular states 
of the hydrostatic theory, whereas panel (c) shows a broken- 
symmetry (phase separated) state. 



from a Maxwellian distribution at temperature To = 1. 
In all simulations the coefficient of normal restitution 
/i = 0.92 and the interior radius R- m t/d = 22.0 were fixed, 
whereas the the number of particles 527 < N < 7800 and 
the aspect ratio 1.5 < f2 < 6 were varied. In terms of 
the three scaled hydrodynamic parameters the heat loss 
parameter A = 81.09 was fixed whereas / and f2 varied. 

To compare the simulation results with predictions 
of our hydrostatic theory, all the measurements were 
performed once the system reached a steady state. This 
was monitored by the evolution of the total kinetic 
energy (1/2) X^=i i which first decays and then, on 
the average, stays constant. 

Steady States. Typical steady state snapshots of the 
system, observed in our MD simulation, are displayed 
in Fig. Panel (a) shows a dilute state where the ra- 
dial density inhomogeneity, though actually present, is 
not visible by naked eye. Panels (b) and (c) do exhibit 
a pronounced radial density inhomogeneity. Apart from 
visible density fluctuations, panels (a) and (b) correspond 
to annular states. Panel (c) depicts a broken-symmetry 
(phase separated) state. When an annular state is ob- 
served, its density profile agrees well with the solution of 
the hydrostatic equations ©-([H]). A typical example of 
such a comparison is shown in Fig. [TJ 

Let us fix the aspect ratio Q of the annulus at not 
too a small value and vary the number of particles N. 
First, what happens on a qualitative level? The simula- 
tions show that, at small N, dilute annular states, similar 
to snapshot (a) in Fig. [51 are observed. As N increases, 
broken-symmetric states start to appear. Well within the 
unstable region, found from hydrodynamics, a high den- 
sity cluster appears, like the one shown in Fig. [5};, and 
performs an erratic motion along the exterior wall. As N 
is increased still further, well beyond the high-/ branch 
of the unstable region, an annular state reappears, as in 
Fig. [5Jd. This time, however, the annular state is denser, 
while its local structure varies from a solid-like (with im- 
perfections such as voids and line defects) to a liquid-like. 

To characterize the spatio-temporal behavior of the 
granulate at a steady state, we followed the position of 
the center of mass (COM) of the granulate. Several ex- 
amples of the COM trajectories are displayed in Fig. [5] 
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FIG. 6: Typical steady state snapshots (the left column) and 
the temporal evolution of the COM (the middle column) and 
of the squared amplitude of the fundamental Fourier mode 
(the right column). The temporal data are sampled each 150 
collisions per particle. Each row corresponds to one simu- 
lation with the indicated parameters. The vertical scale of 
panels a and b was stretched for clarity. 



Here cases (a) and (b) correspond, in the hydrodynamic 
language, to annular states. There are, however, signif- 
icant fluctuations of the COM around the center of the 
annulus. These fluctuations are of course not accounted 
for by hydrodynamic theory. In case (b), where the dense 
cluster develops, the fluctuations are much weaker that 
in case (a). More interesting are cases (c) and (d). They 
correspond to broken-symmetry states: well within the 
phase separation region of the parameter space (case 
c) and close to the phase separation border (case d). 
The COM trajectory in case (c) shows that the gran- 
ular "droplet" performs random motion in the azimuthal 
direction, staying close to the exterior wall. This is in 
contrast with case (d) , where fluctuations are strong both 
in the azimuthal and in the radial directions. Following 
the actual snapshots of the simulation, one observes here 
a very complicated motion of the "droplet" , as well as its 
dissolution into more "droplets" , mergers of the droplets 
etc. Therefore, as in the case of granular phase sepa- 
ration in rectangular geometry [is| . the granular phase 
separation in annular geometry is accompanied by con- 
siderable spatio-temporal fluctuations. In this situation 
a clear distinction between a phase-separated state and 
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FIG. 7: The normalized probability distribution functions 
Pi (^4i/ao) f° r = 3 (the left column) and Q = 5 (the right 
column) for different numbers of particles. 



an annular state, and a comparison between the MD sim- 
ulations and hydrodynamic theory, demand proper diag- 
nostics. We found that such diagnostics are provided by 
the azimuthal spectrum of the particle density and its 
probability distribution. 

Azimuthal Density Spectrum. Let us consider the 
(time-dependent) rescaled density field v(r : 9, t) = 
n(r, 9, t) jn c (where r is rescaled to the interior wall radius 
as before), and introduce the integrated field C>(9,t): 



P(6,t) = J v(r,e,t)rdr. 



(13) 



In a system of N particles, 0(9, t) is normalized so that 

TV 



Tic Rf 



(14) 



c 1 4nt 



Because of the periodicity in 9 the function 0(9, t) can be 



TABLE I: The averaged squared relative amplitudes 
{A\(t))/a%, for the first three modes k = 1, 2 and 3. (a) 



N 


= 2634, n = 


3, (b) N = 5267, = 4, (c) 


N = 1000, 


n 


= 2.25, and (d) N = 1 250, fl 


= 3. 




k 


(a) 


(b) 


(c) 


(d) 


1 


0.66 ± 0.05 


0.39 ±0.04 


0.30 ±0.08 


0.77 ±0.05 


2 


0.04 ± 0.02 


0.05 ± 0.02 


0.07 ±0.01 


0.28 ±0.09 


3 


0.03 ± 0.02 


0.03 ± 0.03 


0.02 ±0.02 


0.11 ±0.08 



expanded in a Fourier series: 

oo 

v{6, t)=a a + Y^ [ofc(t) cos(fc6>) ± b k (t) sm(k9)} , (15) 

k=l 

where ag is independent of time because of the normal- 
ization condition p4[) . We will work with the quantities 

Al(t) = a 2 k (t) + bl(t), k>l. (16) 

For the (deterministic) annular state one has A/~ = for 
all k > 1, while for a symmetry-broken state Ak > 0. 
The relative quantities A^,(t)/a^ can serve as measures 
of the azimuthal symmetry breaking. As is shown in 
Table HI A\ (t) is usually much larger (on the average) 
that the rest of A\{t). Therefore, the quantity A\(t)/dQ 
is sufficient for our purposes. 

Once the system relaxed to a steady state, we followed 
the temporal evolution of the quantity A\/a% ) . Typical 
results are shown in the right column of Fig. [5] One 
observes that, for annular states, this quantity is usu- 
ally small, as is the cases (a) and (b) in Fig. [|5J For 
broken-symmetry states A\ is larger, and it increases as 
one moves deeper into the phase separation region. (No- 
tice that in Fig. [5] the averaged value of A\fa\ in (c) 
is larger than in (d), which means that (c) is deeper in 
the phase separation region.) Another characteristics of 
A\(t)/a^ is the magnitude of fluctuations. One can no- 
tice that, in the vicinity of the phase separation border 
the fluctuations are stronger (as in case (d) in Fig. [S]). 

All these properties are encoded in the probability dis- 
tribution Pi of the values of (Ai/oq) 2 : the ultimate tool 
of our diagnostics. Figure [7] shows two series of measure- 
ments of this quantity at different N: for Q — 3 and 
fi = 5. By following the position of the maximum of 
Pi we were able to to sharply discriminate between the 
annular states and phase separated states and therefore 
to locate the phase separation border. When the max- 
imum of Pi occurs at the zero value of (Ai/oq) 2 (as in 
cases (a) and (d) and, respectively, (e) and (h) in Fig. [7]), 
an annular state is observed. On the contrary, when the 
maximum of Pi occurs at a non-zero value of (Ai/oq) 2 
(as in cases (b) and (c) and, respectively, (f) and (g) 
in Fig. [7]), a phase separated state is observed. In each 
case, the width of the probability distribution (measured, 
for example, at the half-maximum) yields a direct mea- 
sure of the magnitude of fluctuations. Near the phase 







7 / 


1 1 

unstable 


1 




6 




00 ,i / • 


• • 




• 














G4 




d / 










CO •••• • • f 


»• • • • 





- 


2 




!• • 

[m • 

00' 




1 


1 


1 





0.5 


1 


1.5 


2 



f (£2 2 -l)/2 

FIG. 8: The U-f phase diagram for A = 81.09. The solid 
curve is given by the granular hydrostatics: it shows the bor- 
ders of the region where the annular state is unstable with 
respect to small perturbations. The filled symbols depict 
the parameters in which phase separated states are observed, 
whereas the hollow symbols show the parameters at which 
annular states are observed. The dashed line is an estimated 
bin< 




f (£2 2 -l)/2 

FIG. 9: The marginal stability lines for our main setting (the 
solid line) and for an alternative setting in which the thermal 
wall is at r = R ex t and the elastic wall is at r = Ri n t (the 
dashed line). 

separation border, strong fluctuations (that is, broader 
distributions) are observed, as in case (c) of Fig. \7\ 

Using the position of the maximum of Pi as a crite- 
rion for phase separation, we show, in Fig. the ft — f 
diagram obtained from the MD simulations. The same 
figure also depicts the hydrostatic prediction of the in- 
stability region. One can see that the instability region is 
located within the phase separation region, as expected. 

IV. SOME MODIFICATIONS OF THE MODEL 

We also investigated an alternative setting in which the 
exterior wall is the driving wall, while the interior wall 
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Snapshots P } vs. (Aj/aJ 




FIG. 10: Typical steady state snapshots (the left column) and 
the normalized probability distribution functions Pi (Ai/Oq) 
for an inelastic exterior wall, (x wa ll = 0.99 (the right column) 
for different numbers of particles. 



in Fig. [TO] It can be seen that, for the right choice of 
parameters, the phase separation persists. This result is 
important for a possible experimental test of our theory. 
V. SUMMARY 

We combined equations of granular hydrostatics and 
event-driven MD simulations to investigate spontaneous 
phase separation of a monodisperse gas of inelastically 
colliding hard disks in a two-dimensional annulus, the in- 
ner circle of which serves as a "thermal wall" . A marginal 
stability analysis yields a region of the parameter space 
where the annular state - the basic, azimuthally sym- 
metric steady state of the system - is unstable with re- 
spect to small perturbations which break the azimuthal 
symmetry. The physics behind the instability is negative 
effective compressibility of the gas in the azimuthal di- 
rection, which results from the inelastic energy loss. MD 
simulations of this system show phase separation, but 
it is masked by large spatio-temporal fluctuations. By 
measuring the probability distribution of the amplitude 
of the fundamental Fourier mode of the azimuthal spec- 
trum of the particle density we have been able to clearly 
identify the transition to phase separated states in the 
MD simulations. We have found that the instability re- 
gion of the parameter space, predicted from hydrostatics, 
is located within the phase separation region observed 
in the MD simulations. This implies the presence of a 
binodal (coexistence) region, where the annular state is 
metastable, similar to what was found in rectangular ge- 
ometry [lj, [l9[ . The instability persists in an alternative 
setting (a driving exterior wall and an elastic interior 
wall), and also when the elastic wall is replaced by a 
weakly inelastic one. We hope our results will stimulate 
experimental work on the phase separation instability. 



is elastic. The corresponding hydrostatic problem is de- 
termined by the same three scaled parameters /, A and 
fi, but the boundary conditions must be changed accord- 
ingly. Here azimuthally symmetric clusters appear near 
the (elastic) interior wall. Symmetry breaking instabil- 
ity occurs here as well. We found very similar marginal 
stability curves here, but they are narrower (as shown in 
Fig. [5]) than those obtained for the original setting. 

Finally, we returned to our original setting and per- 
formed several MD simulations, replacing the perfectly 
elastic exterior wall by a weakly inelastic one. The in- 
elastic particle collisions with the exterior wall were mod- 
eled in the same way as the inelastic collisions between 
particles. Typical results of these simulations are shown 
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